Title

Data

acc1_cancer_cells = readRDS("./Data/acc1_cancer_cells_15KnCount_V3.RDS")
all_acc_cancer_cells = readRDS("./Data/acc_cancer_cells_V3.RDS")
acc_all_cells = readRDS("./Data/acc_tpm_nCount_mito_no146_15k_with_ACC1_.RDS")
luminal_pathways = c("CHARAFE_BREAST_CANCER_LUMINAL_VS_BASAL_DN","CHARAFE_BREAST_CANCER_LUMINAL_VS_BASAL_UP","CHARAFE_BREAST_CANCER_LUMINAL_VS_MESENCHYMAL_DN","CHARAFE_BREAST_CANCER_LUMINAL_VS_MESENCHYMAL_UP","HUPER_BREAST_BASAL_VS_LUMINAL_DN","LIM_MAMMARY_LUMINAL_PROGENITOR_UP","SMID_BREAST_CANCER_LUMINAL_B_UP" )

# add luminal pathways
luminal_gs = msigdbr(species = "Homo sapiens") %>%as.data.frame() %>% dplyr::filter(gs_name %in% luminal_pathways)%>% dplyr::distinct(gs_name, gene_symbol) %>% as.data.frame()

Parameters

suffix = "01_25"
data_to_read = ""

functions

source_from_github(repositoy = "DEG_functions",version = "0.2.24")

ℹ SHA-1 hash of file is a183df1c565702ecd8ed338bb2abfb0e13415d8e

source_from_github(repositoy = "HMSC_functions",version = "0.1.12",script_name = "functions.R")

ℹ SHA-1 hash of file is 2934dee5f6b9fee69635192f19b1bcc205e05b62

source_from_github(repositoy = "cNMF_functions",version = "0.3.53",script_name = "cnmf_function_Harmony.R")

ℹ SHA-1 hash of file is a6f56985fe95e4026c7830d7100dd81dd950df3b

UMAP

DimPlot(object = acc1_cancer_cells,pt.size = 2)

features UMAP

FeaturePlot(object = acc1_cancer_cells,features = c("TP63","ACTA2","IL12B","CNN1"))

Warning in FeaturePlot(object = acc1_cancer_cells, features = c("TP63",  :
  All cells have the same value (0) of IL12B.

Enrichment analysis HMSC vs ACC

patient.ident = all_acc_cancer_cells$patient.ident %>% as.data.frame()
patient.ident[,1] = as.character(patient.ident[,1])
patient.ident[patient.ident[,1] == "ACC1",] = "HMSC"
patient.ident[,1] = as.factor(patient.ident[,1])
all_acc_cancer_cells = AddMetaData(object = all_acc_cancer_cells,metadata = patient.ident,col.name = "patient.ident")
all_acc_cancer_cells = SetIdent(all_acc_cancer_cells, value ="patient.ident")
acc_deg <- FindMarkers(all_acc_cancer_cells, ident.1 = "HMSC",logfc.threshold = 1.5,features = VariableFeatures(all_acc_cancer_cells))
enrichment_analysis(acc_deg,background = VariableFeatures(all_acc_cancer_cells),fdr_Cutoff = 0.01,ident.1 = "HMSC",ident.2 = "ACC",show_by = 1)

cell cycle filtering

hallmark_name = "GO_MITOTIC_CELL_CYCLE"
genesets  =getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
var_features=all_acc_cancer_cells@assays$RNA@var.features
geneIds= genesets[[hallmark_name]]@geneIds
score <- apply(all_acc_cancer_cells@assays$RNA@scale.data[intersect(geneIds,var_features),],2,mean)
all_acc_cancer_cells=AddMetaData(all_acc_cancer_cells,score,hallmark_name)

#filter:
all_acc_cancer_cells_ccFiltered=all_acc_cancer_cells[,all_acc_cancer_cells@meta.data[[hallmark_name]]< 0.3]


min_threshold = min(all_acc_cancer_cells$GO_MITOTIC_CELL_CYCLE)
max_threshold = max(all_acc_cancer_cells$GO_MITOTIC_CELL_CYCLE)

Before cc filtering

FeaturePlot(object = all_acc_cancer_cells,features = hallmark_name) + ggtitle("Before cc filtering")  & scale_color_gradientn(colours = plasma(n = 10, direction = -1), limits = c(min_threshold, max_threshold))

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

After cc filtering

FeaturePlot(object = all_acc_cancer_cells_ccFiltered,features = hallmark_name) + ggtitle("After cc filtering") & scale_color_gradientn(colours = plasma(n = 10, direction = -1), limits = c(min_threshold, max_threshold))

Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.

Enrichment analysis filtered HMSC vs ACC

patient.ident = all_acc_cancer_cells_ccFiltered$patient.ident %>% as.data.frame()
patient.ident[,1] = as.character(patient.ident[,1])
patient.ident[patient.ident[,1] == "ACC1",] = "HMSC"
patient.ident[,1] = as.factor(patient.ident[,1])
all_acc_cancer_cells_ccFiltered = AddMetaData(object = all_acc_cancer_cells_ccFiltered,metadata = patient.ident,col.name = "patient.ident")
all_acc_cancer_cells_ccFiltered = SetIdent(all_acc_cancer_cells_ccFiltered, value ="patient.ident")
acc_deg <- FindMarkers(all_acc_cancer_cells_ccFiltered, ident.1 = "HMSC",logfc.threshold = 1.5,features = VariableFeatures(all_acc_cancer_cells_ccFiltered))
enrichment_analysis(acc_deg,background = VariableFeatures(all_acc_cancer_cells_ccFiltered),fdr_Cutoff = 0.01,ident.1 = "HMSC",ident.2 = "ACC",show_by = 1)

MYB expression

all_acc_cancer_cells = SetIdent(object = all_acc_cancer_cells,value = "patient.ident") #active snn graph
FeaturePlot(object = all_acc_cancer_cells,features = "MYB",label = T)

CNV

new.cluster.ids <- c("cancer", #0
                     "cancer", #1
                     "CAF", #2
                     "cancer", #3
                     "Endothelial", #4
                     "cancer", #5
                     "cancer", #6
                     "CAF", #7
                     "CAF", #8
                     "CAF", #9
                     "cancer", #10
                     "CAF", #11
                     "cancer", #12
                     "cancer", #13
                     "cancer", #14
                     "cancer", #15
                     "cancer", #16
                     "WBC", #17
                     "CAF" #18
                     )

#rename idents:
acc_all_cells = SetIdent(object = acc_all_cells,value = "RNA_snn_res.1") #active snn graph
names(new.cluster.ids) <- levels(acc_all_cells) #add snn graph levels to new.cluster.ids
acc_all_cells@meta.data[["seurat_clusters"]] = acc_all_cells@meta.data[["RNA_snn_res.1"]]
acc_all_cells = SetIdent(object = acc_all_cells,value = "seurat_clusters")
acc_all_cells <- RenameIdents(acc_all_cells, new.cluster.ids) 

# divide "cancer" into patients:
cell_types = acc_all_cells@active.ident %>% as.data.frame()
cell_types[,1]<- as.character(cell_types[,1])
cell_types = cbind(cell_types,acc_all_cells$patient.ident) %>% setnames(old = names(.), 
         new = c('cell_type','patient'))
cell_types[cell_types$cell_type == "cancer",] = cell_types[cell_types$cell_type == "cancer",2]


# hmsc_rows = (startsWith(x = rownames(cell_types),prefix = "ACC.plate2") | startsWith(x = rownames(cell_types),prefix = "ACC1.")) & cell_types[,1] == "cancer" 
# acc_rows = !(startsWith(x = rownames(cell_types),prefix = "ACC.plate2") | startsWith(x = rownames(cell_types),prefix = "ACC1.")) & cell_types[,1] == "cancer" 
# cell_types[,1][hmsc_rows]  = "HMSC"
# cell_types[,1][acc_rows]  = "ACC"

#add to metadata:
cell_types[,2] = NULL 
cell_types[cell_types$cell_type == "ACC1",] = "HMSC"
acc_all_cells = AddMetaData(object =acc_all_cells ,metadata = cell_types,col.name = "cell.type")

CNV UMAP

CNV plot

## {-}

CNV subtypes

cnv_subtypes = as.data.frame(cutree(infercnv_obj_default@tumor_subclusters[["hc"]][["HMSC"]], k = 2))
names(cnv_subtypes)[1] = "cnv.cluster"
rownames(cnv_subtypes) = rownames(cnv_subtypes) %>% gsub(pattern = "-",replacement = "\\.")
infercnv.observations = data.frame(fread(file = "./Data/inferCNV/infercnv.observations.txt"), row.names=1)

Warning in fread(file = "./Data/inferCNV/infercnv.observations.txt") :
  Detected 1332 column names but the data has 1333 columns (i.e. invalid file). Added 1 extra default column name for the first column which is guessed to be row names or an index. Use setnames() afterwards if this guess is not correct, or fix the file write command that created the file to create a valid file.

names_to_keep = colnames(infercnv.observations) %in% (colnames(acc1_cancer_cells) %>% gsub(pattern = "_",replacement = "\\."))
infercnv.observations = infercnv.observations[,names_to_keep]

rotate <- function(x) t(apply(x, 2, rev))
infercnv.observations2 = infercnv.observations %>% rotate() %>%  rotate() %>% rotate()%>% as.data.frame() 
breaks = c(0.700891861704857,
0.742366945528369,
0.783842029351881,
0.825317113175393,
0.866792196998905,
0.908267280822417,
0.949742364645928,
0.99121744846944,
1.03269253229295,
1.07416761611646,
1.11564269993998,
1.15711778376349,
1.198592867587,
1.24006795141051,
1.28154303523402,
1.32301811905753)
pheatmap(infercnv.observations2,cluster_cols = F,cluster_rows = F, show_rownames = F,show_colnames = F, breaks = breaks,color = colorRampPalette(rev(c("darkred", "white", "darkblue")))(15),annotation_row = cnv_subtypes)

rownames(cnv_subtypes) = rownames(cnv_subtypes) %>% gsub(pattern = "2\\.",replacement = "2_")
rownames(cnv_subtypes) = rownames(cnv_subtypes) %>% gsub(pattern = "3\\.",replacement = "3_")
acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = cnv_subtypes)

DimPlot(acc1_cancer_cells,group.by = "cnv.cluster",pt.size = 2,cols =colors)

Original score

original_myo_genes = c( "TP63", "TP73", "CAV1", "CDH3", "KRT5", "KRT14", "ACTA2", "TAGLN", "MYLK", "DKK3")
original_lum_genes = c("KIT", "EHF", "ELF5", "KRT7", "CLDN3", "CLDN4", "CD24", "LGALS3", "LCN2", "SLPI" )

calculate_score(dataset = all_acc_cancer_cells,myo_genes = original_myo_genes,lum_genes = original_lum_genes)

correlation of lum score and myo score: -0.45
correlation of lum score and original lum score: 1
correlation of myo score and original myo score: 1

Original score of ACC1

calculate_score(dataset = acc1_cancer_cells,myo_genes = original_myo_genes,lum_genes = original_lum_genes,lum_threshold = 0,myo_threshold = 0)

correlation of lum score and myo score: 0.06
correlation of lum score and original lum score: 1
correlation of myo score and original myo score: 1

0.35 Most correlated score

Myo genes

myo_protein_markers = c("CNN1", "TP63","ACTA2")
top_myo  = top_correlated(dataset = acc1_cancer_cells, genes = myo_protein_markers,threshold = 0.35)
print("Number of genes = " %>% paste(length(top_myo)))

[1] "Number of genes =  20"

message("Names of genes:")

Names of genes:

top_myo %>% head(30)

 [1] "COL16A1"     "RP1-39G22.4" "ACTG2"       "CD200"       "MYLK"        "TP63"        "KCNMB1"     
 [8] "ADAMTS2"     "LOXL2"       "TPM2"        "CLIC3"       "SNCG"        "ACTA2"       "TAGLN"      
[15] "A2M"         "NGFR"        "CNN1"        "PPP1R14A"    "MYL9"        "POM121L9P"  

message("Genes that also apeared in the original score:")

Genes that also apeared in the original score:

base::intersect(top_myo,original_myo_genes) 

[1] "MYLK"  "TP63"  "ACTA2" "TAGLN"

myo_enrich_res = genes_vec_enrichment(genes = top_myo,background = rownames(acc1_cancer_cells),homer = T,title = "myo top enrichment",custom_pathways = luminal_gs)

myo_enrich_res

Lum genes

lum_protein_markers = c("KIT")
top_lum  = top_correlated(dataset = acc1_cancer_cells, genes = lum_protein_markers,threshold = 0.35,n_vargenes = 5000)

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

print("Number of genes = " %>% paste(length(top_lum)))

[1] "Number of genes =  5"

message("Names of genes:")

Names of genes:

top_lum %>% head(30)

[1] "B3GNT2"  "GLRB"    "EFNA5"   "ALDH3B2" "CCND1"  

message("Genes that also apeared in the original score:")

Genes that also apeared in the original score:

base::intersect(top_lum,original_lum_genes) 

character(0)

lum_enrich_res = genes_vec_enrichment(genes = top_lum,background = rownames(acc1_cancer_cells),homer = T,title = "lum top enrichment",custom_pathways = luminal_gs)

lum_enrich_res

top correlated score

calculate_score(dataset = acc1_cancer_cells,myo_genes = top_myo,lum_genes = top_lum,lum_threshold = 0,myo_threshold = 0)

correlation of lum score and myo score: -0.05
correlation of lum score and original lum score: 0.58
correlation of myo score and original myo score: 0.77

enriched genes score

rownames(lum_enrich_res) = lum_enrich_res$pathway_name
lum_enriched_genes = lum_enrich_res[1,"geneID"] %>% strsplit(split = "/") %>% .[[1]] %>% c(.,lum_protein_markers) #add original markers

rownames(myo_enrich_res) = myo_enrich_res$pathway_name
myo_enriched_genes = myo_enrich_res[1,"geneID"] %>% strsplit(split = "/") %>% .[[1]] %>% c(.,myo_protein_markers) #add original markers

calculate_score(dataset = acc1_cancer_cells,myo_genes = myo_enriched_genes,lum_genes = lum_enriched_genes,lum_threshold = -2.5,myo_threshold = -2.5)

correlation of lum score and myo score: -0.16
correlation of lum score and original lum score: 0.43
correlation of myo score and original myo score: 0.79

0.2 Most correlated score

myo Genes

n_vargenes = 2000
myo_protein_markers = c("CNN1", "TP63","ACTA2")
top_myo  = top_correlated(dataset = acc1_cancer_cells, genes = myo_protein_markers,threshold = 0.2,n_vargenes = n_vargenes)
print("Number of genes = " %>% paste(length(top_myo)))

[1] "Number of genes =  48"

message("Names of genes:")

Names of genes:

top_myo %>% head(30)

 [1] "PLOD1"         "RP1-39G22.4"   "PDZK1"         "CSRP1"         "CHI3L1"        "HIST3H3"      
 [7] "AC104699.1"    "ACTG2"         "LYG1"          "DALRD3"        "CD200"         "RP11-627C21.1"
[13] "FGFBP2"        "IGFBP7"        "HSPB3"         "RP11-168A11.4" "SPARC"         "CD83"         
[19] "HLA-DOB"       "TBCC"          "NFKBIE"        "MIR3662"       "SMOC2"         "FBXL6"        
[25] "TPM2"          "SNCG"          "ACTA2"         "DKK3"          "MIR7113"       "CTA-797E19.3" 

message("Genes that also apeared in the original score:")

Genes that also apeared in the original score:

base::intersect(top_myo,original_myo_genes) 

[1] "ACTA2" "DKK3"  "TAGLN"

myo_enrich_res = genes_vec_enrichment(genes = top_myo,background = VariableFeatures(acc1_cancer_cells) %>% head(n_vargenes),homer = T,title = "myo top enrichment",custom_pathways = luminal_gs)

myo_enrich_res

Lum Genes

acc1_cancer_cells = FindVariableFeatures(object = acc1_cancer_cells,nfeatures = 15000)

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

lum_protein_markers = c("KIT")
n_vargenes = 2000
top_lum  = top_correlated(dataset = acc1_cancer_cells, genes = lum_protein_markers,threshold = 0.3,n_vargenes = n_vargenes)
print("Number of genes = " %>% paste(length(top_lum)))

[1] "Number of genes =  13"

message("Names of genes:")

Names of genes:

top_lum %>% head(30)

 [1] "CSDE1"   "RRP9"    "TKT"     "TFG"     "ATP1B3"  "EFNA5"   "GABRP"   "CALML5"  "MMP7"    "SNORD68"
[11] "APP"     "SDF2L1"  "SAT1"   

message("Genes that also apeared in the original score:")

Genes that also apeared in the original score:

base::intersect(top_lum,original_lum_genes) 

character(0)

lum_enrich_res = genes_vec_enrichment(genes = top_lum,background = VariableFeatures(acc1_cancer_cells) %>% head(n_vargenes),homer = T,title = "lum top enrichment",custom_pathways = luminal_gs)

lum_enrich_res

calculate_score(dataset = acc1_cancer_cells,myo_genes = top_myo,lum_genes = top_lum,lum_threshold = 2,myo_threshold = 1)

correlation of lum score and myo score: -0.02
correlation of lum score and original lum score: 0.65
correlation of myo score and original myo score: 0.77

myo_intersected = intersect(top_myo,original_myo_genes) 
lum_intersected = intersect(top_lum,original_lum_genes) 
message("genes in myo score:")

genes in myo score:

myo_intersected

[1] "ACTA2" "DKK3"  "TAGLN"

message("genes in lum score:")

genes in lum score:

lum_intersected

[1] "LGALS3"

calculate_score(dataset = acc1_cancer_cells,myo_genes = myo_intersected,lum_genes = lum_intersected,lum_threshold = 2,myo_threshold = 1)

correlation of lum score and myo score: 0.03
correlation of lum score and original lum score: 0.69
correlation of myo score and original myo score: 0.81

enriched genes

rownames(lum_enrich_res) = lum_enrich_res$pathway_name
lum_enriched_genes = lum_enrich_res[3,"geneID"] %>% strsplit(split = "/") %>% .[[1]] %>% c(.,lum_protein_markers) #add original markers

rownames(myo_enrich_res) = myo_enrich_res$pathway_name
myo_enriched_genes = myo_enrich_res[3,"geneID"] %>% strsplit(split = "/") %>% .[[1]] %>% c(.,myo_protein_markers) #add original markers

message("genes in myo score:")

genes in myo score:

myo_enriched_genes

 [1] "FBLIM1"  "NEXN"    "NMNAT2"  "MYLK"    "CCDC50"  "IGFBP7"  "RAI14"   "ARAP3"   "SPARC"   "CALD1"   "LOXL2"  
[12] "COL5A1"  "ACTA2"   "DKK3"    "MSRB3"   "COL4A1"  "ACTN1"   "TPM1"    "TGFB1I1" "ADORA2B" "MXRA7"   "CNN1"   
[23] "TP63"    "ACTA2"  

message("genes in lum score:")

genes in lum score:

lum_enriched_genes

 [1] "PADI2"      "PATJ"       "PEX11B"     "APH1A"      "C1orf43"    "EFNA4"      "NECTIN4"    "SCYL3"      "ELF3"      
[10] "SOX13"      "IRF6"       "MED28"      "EPB41L4A"   "RGL2"       "C6orf132"   "TPD52L1"    "ICA1"       "MACC1"     
[19] "TRPS1"      "FAM83H"     "RASEF"      "ARRDC1"     "COMMD3"     "ANK3"       "GSTO2"      "PDCD4"      "EHF"       
[28] "ALDH3B2"    "SHANK2"     "SORL1"      "FKBP4"      "PTPN6"      "DUSP16"     "RERG"       "ADCY6"      "ERBB3"     
[37] "ERP29"      "SUSD6"      "RPS6KA5"    "SPINT1"     "FEM1B"      "TLE3"       "SCAMP2"     "CLN3"       "ADGRG1"    
[46] "ATP2C2"     "GGT6"       "MYO1D"      "ST6GALNAC2" "CYB5A"      "BLVRB"      "VRK3"       "SYCP2"      "TMPRSS2"   
[55] "LIMK2"      "KIT"       

calculate_score(dataset = acc1_cancer_cells,myo_genes = myo_enriched_genes,lum_genes = lum_enriched_genes,lum_threshold = 0,myo_threshold = -1)

correlation of lum score and myo score: 0.1
correlation of lum score and original lum score: 0.71
correlation of myo score and original myo score: 0.76

enriched genes and in original score

myo_enriched_genes = myo_enriched_genes[myo_enriched_genes %in% original_myo_genes]
lum_enriched_genes = lum_enriched_genes[lum_enriched_genes %in% original_lum_genes]

message("genes in myo score:")

genes in myo score:

myo_enriched_genes

[1] "MYLK"  "ACTA2" "DKK3"  "TP63"  "ACTA2"

message("genes in lum score:")

genes in lum score:

lum_enriched_genes

[1] "EHF" "KIT"

calculate_score(dataset = acc1_cancer_cells,myo_genes = myo_enriched_genes,lum_genes = lum_enriched_genes,lum_threshold = 2,myo_threshold = 2)

correlation of lum score and myo score: -0.07
correlation of lum score and original lum score: 0.62
correlation of myo score and original myo score: 0.77

HPV

Only HMSC cancer cells:

HPV33_P3 = fread("./Data/HPV33_P3.txt",col.names = c("plate","reads")) %>% as.data.frame()
HPV33_P3.df = HPV33_P3 %>% mutate(
  plate = gsub(x =HPV33_P3$plate, replacement = "",pattern = "_.*$") 
  %>% gsub(pattern = "-P",replacement = ".P") 
  %>% gsub(pattern = "-",replacement = "_",)
  )
HPV33_P3.df = HPV33_P3.df %>% dplyr::filter(HPV33_P3.df$plate %in% colnames(acc1_cancer_cells))
rownames(HPV33_P3.df)  <- HPV33_P3.df$plate
HPV33_P3.df$plate = NULL


HPV33_P2 = fread("./Data/HPV33_P2.txt",col.names = c("plate","reads")) %>% as.data.frame()
HPV33_P2.df = HPV33_P2 %>% mutate(
  plate = gsub(x =HPV33_P2$plate, replacement = "",pattern = "_.*$") 
  %>% gsub(pattern = "plate2-",replacement = "plate2_",)
  %>% gsub(pattern = "-",replacement = "\\.",)
  )
HPV33_P2.df = HPV33_P2.df %>% dplyr::filter(HPV33_P2.df$plate %in% colnames(acc1_cancer_cells))
rownames(HPV33_P2.df)  <- HPV33_P2.df$plate
HPV33_P2.df$plate = NULL

HPV33 = rbind(HPV33_P3.df,HPV33_P2.df)
acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = HPV33,col.name = "HPV33.reads")
FeaturePlot(acc1_cancer_cells,features = "HPV33.reads",max.cutoff = 40)

data = FetchData(object = acc1_cancer_cells,vars = "HPV33.reads")
print(
  data %>% 
  ggplot(aes( x=HPV33.reads)) + 
  geom_density()
)

hpv33_positive = HPV33 %>% dplyr::mutate(hpv33_positive = case_when(reads >= 750 ~"strong positive",
                                                                    reads < 750 & reads > 10 ~ "positive",
                                                                    reads < 10 ~ "negative")
)



hpv33_positive$reads = NULL
acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = hpv33_positive)

DimPlot(object = acc1_cancer_cells,group.by  = c("hpv33_positive"),pt.size = 2)

cNMF

library(reticulate)

#write expression
acc1_cancer_cells = FindVariableFeatures(object = acc1_cancer_cells,nfeatures = 2000)
vargenes = VariableFeatures(object = acc1_cancer_cells)
hmsc_expression = t(as.matrix(GetAssayData(acc1_cancer_cells,slot='data')))
hmsc_expression = 2**hmsc_expression #convert from log2(tpm+1) to tpm
hmsc_expression = hmsc_expression-1
# hmsc_expression = hmsc_expression[,!colSums(hmsc_expression==0, na.rm=TRUE)==nrow(hmsc_expression)] #delete rows that have all 0
hmsc_expression = hmsc_expression[,vargenes]
write.table(x = hmsc_expression ,file = './Data/cNMF/hmsc_expressionData_2Kvargenes.txt',sep = "\t")

from cnmf import cNMF
name = 'HMSC_cNMF_2Kvargenes'
outdir = './Data/cNMF'
K_range = np.arange(3,10)
cnmf_obj = cNMF(output_dir=outdir, name=name)
counts_fn='./Data/cNMF/hmsc_expressionData_2Kvargenes.txt'
tpm_fn = counts_fn ## This is a weird case where because this dataset is not 3' end umi sequencing, we opted to use the TPM matrix as the input matrix rather than the count matrix

cnmf_obj.prepare(counts_fn=counts_fn, components=K_range, seed=14,tpm_fn=tpm_fn)

cnmf_obj.factorize(worker_i=0, total_workers=1)

cnmf_obj.combine()
cnmf_obj.k_selection_plot()

Save object

import pickle
f = open('./Data/cNMF/HMSC_cNMF_2Kvargenes/cnmf_obj.pckl', 'wb')
pickle.dump(cnmf_obj, f)
f.close()

Load object

from cnmf import cNMF
import pickle
f = open('./Data/cNMF/HMSC_cNMF_2Kvargenes/cnmf_obj.pckl', 'rb')
cnmf_obj = pickle.load(f)
f.close()

selected_k = 4
density_threshold = 0.1
cnmf_obj.consensus(k=selected_k, density_threshold=density_threshold,show_clustering=True)

/sci/labs/yotamd/lab_share/avishai.wizel/python_envs/miniconda/envs/cnmf_env_6/lib/python3.7/site-packages/scanpy/preprocessing/_simple.py:843: UserWarning: Received a view of an AnnData. Making a copy.
  view_to_actual(adata)

usage_norm, gep_scores, gep_tpm, topgenes = cnmf_obj.load_results(K=selected_k, density_threshold=density_threshold)

# calculate usage by expression*genes coefs:
gep_scores = py$gep_scores
all_metagenes= py$usage_norm

Enrichment analysis by top 200 genes of each program

plt_list = list()
for (i in 1:ncol(gep_scores)) {
  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),200) #take top top_genes_num
  res = genes_vec_enrichment(genes = top,background = rownames(gep_scores),homer = T,title = 
                    i,silent = T,return_all = T)
   
  plt_list[[i]] = res$plt
}
gridExtra::grid.arrange(grobs = plt_list)

plt_list = list()
for (i in 1:ncol(gep_scores)) {
  top_genes = gep_scores  %>%  arrange(desc(gep_scores[i])) #sort by score a
  top = head(rownames(top_genes),10) #take top top_genes_num
  message(paste("program ",i,"top genes:"))
  print(top)

}

program  1 top genes:

 [1] "IGKV5-2" "NDUFB7"  "LGALS1"  "UBL5"    "S100A6"  "S100A11" "CUTA"    "IFITM3"  "JTB"     "IL17RB" 

program  2 top genes:

 [1] "EGR1"    "C6orf62" "JUNB"    "DNAJA1"  "ATF3"    "IER2"    "ERRFI1"  "KLF6"    "CDKN1A"  "MTHFD2" 

program  3 top genes:

 [1] "RP1-128M12.3"  "RP11-374F3.3"  "RP11-403A21.3" "RP11-454L9.2"  "AC097500.1"    "RP11-463H24.4"
 [7] "RP11-515O17.2" "RP11-536L3.4"  "PALD1"         "AL049758.2"   

program  4 top genes:

 [1] "LTF"           "PIGR"          "FMO2"          "HSD17B2"       "MLPH"          "PRR15L"       
 [7] "RF00019.219"   "RP11-817J15.2" "CD14"          "CLU"          

# Make metagene names
for (i in 1:ncol(all_metagenes)) {
  colnames(all_metagenes)[i] = "metagene." %>% paste0(i)
}

#add each metagene to metadata
for (i in 1:ncol(all_metagenes)) {
  metage_metadata = all_metagenes %>% select(i)
  acc1_cancer_cells = AddMetaData(object = acc1_cancer_cells,metadata = metage_metadata)
}

Note: Using an external vector in selections is ambiguous.
ℹ Use `all_of(i)` instead of `i` to silence this message.
ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
This message is displayed once per session.

FeaturePlot(object = acc1_cancer_cells,features = colnames(all_metagenes))

assignment UMAP

larger_by = 2
message(paste("larger_by = ", larger_by))

larger_by =  2

acc1_cancer_cells = program_assignment(dataset = acc1_cancer_cells,larger_by = larger_by,program_names = colnames(all_metagenes))
selected_k =4
colors =  rainbow(selected_k)
colors = c(colors,"grey")
DimPlot(acc1_cancer_cells,group.by = "program.assignment",pt.size = 2,cols =colors)

show cell cycle program:

hallmark_name = "GO_MITOTIC_CELL_CYCLE"
genesets  =GSEABase::getGmt("./Data/h.all.v7.0.symbols.pluscc.gmt")
var_features=acc1_cancer_cells@assays$RNA@var.features
geneIds= genesets[[hallmark_name]]@geneIds
score <- apply(acc1_cancer_cells@assays$RNA@data[intersect(geneIds,var_features),],2,mean)
acc1_cancer_cells=AddMetaData(acc1_cancer_cells,score,hallmark_name)

FeaturePlot(object = acc1_cancer_cells,features = hallmark_name)

cc_vs_program2 = FetchData(object = acc1_cancer_cells,vars = c("metagene.2",hallmark_name))
cor(cc_vs_program2[1],cc_vs_program2[2])

           GO_MITOTIC_CELL_CYCLE
metagene.2              0.611673

Comparisions

cnv_vs_hpv = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","hpv33_positive"))
test <- fisher.test(table(cnv_vs_hpv))
ggbarstats(
  cnv_vs_hpv, cnv.cluster, hpv33_positive,
  results.subtitle = FALSE,
  subtitle = paste0(
    "Fisher's exact test", ", p-value = ",
    ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
  )
)

cnv_vs_cnmf = FetchData(object = acc1_cancer_cells,vars = c("program.assignment","cnv.cluster"))
cnv_vs_cnmf = cnv_vs_cnmf %>% dplyr::filter(program.assignment == "1" |program.assignment == "2" )
test <- fisher.test(table(cnv_vs_cnmf))
ggbarstats(
  cnv_vs_cnmf, program.assignment, cnv.cluster,
  results.subtitle = FALSE,
  subtitle = paste0(
    "Fisher's exact test", ", p-value = ",
    ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
  )
)

cnmf_vs_hpv = FetchData(object = acc1_cancer_cells,vars = c("program.assignment","hpv33_positive"))
cnmf_vs_hpv = cnmf_vs_hpv %>% dplyr::filter(program.assignment == "1" |program.assignment == "2" )
test <- fisher.test(table(cnmf_vs_hpv))
ggbarstats(
  cnmf_vs_hpv, program.assignment, hpv33_positive,
  results.subtitle = FALSE,
  subtitle = paste0(
    "Fisher's exact test", ", p-value = ",
    ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 3))
  )
)

myb_vs_cnv = FetchData(object = acc1_cancer_cells,vars = c("cnv.cluster","MYB"))
myb_vs_cnv $cnv.cluster = as.character(myb_vs_cnv $cnv.cluster )

ggboxplot(myb_vs_cnv, x = "cnv.cluster", y = "MYB",
          palette = "jco",
          add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("1","2")))

myb_vs_hpv = FetchData(object = acc1_cancer_cells,vars = c("hpv33_positive","MYB"))
myb_vs_hpv $hpv33_positive = as.character(myb_vs_hpv $hpv33_positive )

ggboxplot(myb_vs_hpv, x = "hpv33_positive", y = "MYB",
          palette = "jco",
          add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("positive","negative"),c("strong positive", "positive"),c("strong positive", "negative")))+ stat_summary(fun.data = function(x) data.frame(y=15, label = paste("Mean=",round(mean(x),digits = 2))), geom="text") +ylab("log2(MYB)")

Warning in wilcox.test.default(c(7.30445531374182, 0, 7.65159269158999,  :
  cannot compute exact p-value with ties

hpvReads_vs_myb = FetchData(object = acc1_cancer_cells,vars = c("HPV33.reads","MYB"))
corr = cor(hpvReads_vs_myb$HPV33.reads,hpvReads_vs_myb$MYB)
print("correlation of MYB abd hpv33_reads:" %>% paste(corr %>% round(digits = 2)) )

[1] "correlation of MYB abd hpv33_reads: 0.21"